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Abstract 



Some general dynamical properties of models for compaction of granular me- 
dia based on master equations are analyzed. In particular, a one-dimensional 
lattice model with short-ranged dynamical constraints is considered. The 
stationary state is consistent with Edward's theory of powders. The sys- 
tem is submitted to processes in which the tapping strength is monotonically 
increased and decreased. In such processes the behavior of the model resem- 
bles the reversible-irreversible branches which have been recently observed in 
experiments. This behavior is understood in terms of the general dynamical 
properties of the model, and related to the hysteresis cycles exhibited by struc- 
tural glasses in thermal cycles. The existence of a "normal" solution, i. e., a 
special solution of the master equation which is monotonically approached by 
all the other solutions, plays a fundamental role in the understanding of the 
hysteresis effects. 

PACS numbers: 81.05.Rm, 05.50.+q, 81.20.Ev 

Keywords: granular media, master equation, hysteresis, normal solution. 
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I. INTRODUCTION 



The phenomenology of granular materials is very rich. One of the most characteristic 
dynamical behaviors is compaction. If a system of grains in a loosely packed configuration 
is submitted to vertical vibrations, or tappings, it approaches slowly a more compact state. 
The relaxation is well described by an inverse logarithmic law, and both the relaxation 
function and the density of the "stationary" state seem to depend only on a dimensionless 
parameter characterizing the vibration intensity [|TJ. 

More recently, it has also been shown that granular materials exhibit irreversible- 
reversible cycles if the vibration strength is increased and decreased alternately. When 
the system is tapped with increasing intensity starting from a loosely packed state, the den- 
sity shows a non-monotonic behavior, with a maximum at a certain value of the intensity of 
vibration. If, afterwards, the granular medium is tapped with decreasing intensity, the den- 
sity increases monotonically, and hysteresis effects show up. Interestingly, if the vibration 
intensity is again increased, the density follows a curve that is approximately equal to the 
evolution in the decreasing intensity process, being thus "reversible" p],|3[]. Of course, the 
rate of variation of the vibration intensity is the same for all the processes described above. 

The static properties of granular materials have been studied by Edwards and co-workers 
0|J. The starting point is the plausible idea that all the "microscopic" configurations of 
a powder having the same volume are equiprobable, provided that the powder has been 
prepared by extensive manipulation, i. e., by processes which do not act on individual grains. 
This has been called the "ergodic hypothesis" of powders. On this basis, it is possible to make 
an analogy between the variables of a molecular system and the parameters characterizing 
the state of a powder. The volume of a powder is analogous to the energy, while the entropy 
remains the same quantity, measuring the available number of configurations. The derivative 
of the volume with respect to the entropy is called the "compact ivity" of the powder, and 
plays the role of the temperature. The lowest compactivity corresponds to the densest 
state (minimum volume), and the highest compactivity to the fluffiest stable configuration 
(maximum volume). The stationary value of the volume is an increasing function of the 
compactivity, like the energy increases with the temperature in a thermal system. 

It seems interesting to try to establish a connection between the dynamical and the 
equilibrium properties of powders. It looks reasonable that, if the steady state is reached 
in a tapping process, there should be a relationship between the vibration intensity and 
the compactivity of a powder. The tapping process allows the system to explore the phase 
space of available configurations, and the same role is played by the temperature in a thermal 
system. In this way, it is tempting to explain the tendency of a granular system to move 
over an almost reversible curve, when the tapping strength varies in time, as the approach 
towards the stationary curve of the powder. If the above is true, the compactivity will be 
an increasing function of the vibration intensity, because the density over the "reversible" 
curve decreases with the tapping strength. In fact, this kind of behavior has been found in 
a simple model of a granular system described in terms of a Fokker-Planck equation M. 

Also, glass formers exhibit hysteresis effects when cooled and heated through the glass 
transition region. When a glass former is cooled, the system follows the equilibrium curve 
until a certain temperature T g , at which it gets frozen due to the fast increase of the relax- 
ation time. If this structural glass is reheated from its frozen state, the equilibrium curve 
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is only approached for temperatures larger than T g . This is due to the fact that the system 
starts from a configuration in which the structural rearrangements are very difficult [0]. In 
this way, also hysteresis effects show up in glass formers, when they are submitted to thermal 
cycles. 

The study of simple models have been very useful in order to understand, at least qual- 
itatively, the behavior of structural glasses. In particular, hysteresis cycles are also shown 
by simple models when cooled and heated The analytical approach to thermal cycles 

is a very difficult task, because it is necessary to solve the kinetic equations of the model 
with time-dependent coefficients. Nevertheless, some models whose dynamics is formulated 
by means of a simple master equation can be exactly solved ||10|,|rTJ. In those situations, 
the role played by the "normal solution" of the master equation, which is monotonically 



approached by all the other solutions [|12[], has been shown to play a fundamental role. 

Here we will consider a simple model for granular media, which has been previously in- 
troduced ||13||. Its dynamical evolution is governed by a master equation, with the transition 



rates given as functions of the tapping strength. When the system is vibrated with a given 
intensity from a low density state, the density of the system increases monotonically until 
reaching an stationary value. The relaxation of the density is very slow, being very well fit- 
ted by an empirical inverse-logarithmic law. Interestingly, the stationary state of the system 
is described by Edward's theory of powders. Then, it seems worth studying its behavior 
when the tapping strength changes in time for several reasons. Firstly, in order to verify if 
its behavior resembles that of real granular systems, showing the irreversible and reversible 
branches found in experiments 0. Secondly, to understand whether the hysteresis effects are 
related to the existence of a "normal solution" of the master equation with time- dependent 
tapping strength. 

The paper is organized as follows. In Sec. [TI| some general properties of models for 
granular media based on a master equation formulation of the dynamics are considered. 
The existence of a "normal solution" of the master equation is analyzed for the case of 
time-dependent transition rates. Also, the conditions to be verified by the law of variation 
of the tapping strength in order to guarantee that the system approaches the steady state 
curve are discussed. Section [TT1] is devoted to the analysis of the normal solution by means 
of Hilbert's expansion of the master equation. A quite general expression for the first order 
deviation of the normal solution from the stationary curve is obtained. The specific model 
to be considered is presented is Sec. [TV], as well as a brief discussion of its equilibrium 
state. Section [VJ deals with the behavior of the model in processes with time-dependent 
vibration intensity. Firstly, processes in which the tapping strength decreases in time are 
studied. The existence of a phenomenon similar to the laboratory glass transition is analyzed. 
Secondly, we discuss processes with increasing vibration intensity. The role played by the 
normal solution turns out to be fundamental. Finally, the main conclusions of the paper are 
summarized in Sec. [VI]. 

II. SOME GENERAL DYNAMICAL PROPERTIES 

Let us consider a model system whose dynamics is described by means of the master 
equation W4\ 
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dpi(t) 
dt 



(i) 



Here Pi(t) is the probability of finding the system in state % at time t, and the rates Wij for 
transition from state j to state i depend on time in a given way, independently of the state 
of the system. Let us define a function 



Pitt) 
P'i(t) 



(2) 



where Pi(t) and p[(t) are two solutions of Eq. (fll) corresponding to different initial conditions. 
The above definition for H(t) assumes that p'^t) is positive for all the states i. This condition 
will be fulfilled after a transient period if the process defined by Eq. ([]]) is irreducible, even 
in the case that some initial probabilities vanish. The time variation of H(t) is given by 



dH(t) 
dt 



A(t) 



with A(t) being a complicated functional of the two solutions p and p' 
written in the form 



(3) 



that can be 



A{t) = Y,W l3 p' j 



J>j Pi) \ Pi 



, Pi, Pi Pj, Pj 

+ - In - - -i In -i 

Pi Pi Pj Pj 



(4) 



Its main property is that A(t) < 0. Besides, if the transition rates define an irreducible 
process at time t, the equality sign holds only when 



Pi{t)=Pi(t) 



(5) 



for all the states i. As the function H(t) is bounded below, H(t) > 0, it must tend to a 
limit and, therefore, 



lim Pi(t) = lim 



(6) 



Thus all the solutions of the master equation converge toward the same behavior, if the long- 
time limit of the transition rates still define an irreducible process. This equation can be 
understood as showing the existence of a long-time regime, where the influence of the initial 
conditions has been lost, and the state of the system and its dynamics is fully determined by 
the law of variation of the transition rates. Therefore, there will be a special solution of the 
master equation such that all the other solutions approach it after an initial transient period. 
We will refer to this special solution as the "normal" solution of the master equation for the 
given time dependence of the transition rates. A more detailed and general discussion of the 
H-theorem leading to the existence of this 
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'normal" solution can be found in Ref. 
Now suppose that the master equation models the dynamical behavior of a granular 
system submitted to vertical vibration. The transition rates would be functions of the 
parameter T characterizing the strength of the vibration. If the granular pile is vibrated 
with sinusoidal pulses of amplitude A and frequency ui, the quantity T = Auj 2 / g 
where g is the gravity, is usually defined. In the case of time dependent intensity T, the 
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equation determining the evolution of the granular pile will be of the type given by Eq. ([[]). 
Assuming that for arbitrary r ^ the tapping process allows the system to explore the 
whole configuration space, the stochastic process will be irreducible and the existence of a 
normal solution for a given program of variation of the intensity T follows, provided that V 
does not vanish in the long-time limit. 

Let us assume that for every value of the intensity T > 0, the equation 



£^(r>(. s) (r) = E^( r k (sj (r) 



(7) 



has a "canonical" solution, of the form 



rf°(r) 



Z{X) 



exp 



V, 



where Vi is the volume of the system in state i, 

Z(X) = £exp 



AX 



AX 



(8) 



(9) 



is a partition function, A is a constant with the dimension of volume, and X is a variable 
termed the compactivity of the granular system in the framework of Edward's statistical 
mechanics theory of powders [|]||. Of course, in the context considered here the compactivity 
X will be a function of the vibration intensity T, X = f(T), which has to be found for each 
particular model. Equation (^) defines the steady state reached by the system if the vibration 
intensity is constant in time. The macroscopic value of the volume in this state is 



v (s \x) 



(10) 



and the configurational entropy S(X) reads 

S(X) = -\Y / p ( f\x)lnp^(X) 



(11) 



Of course, the compactivity X can be obtained from its usual definition, 



X 



dV (s) 
dS 



(12) 



For the sake of simplicity, in the following we will take the unit of volume such that A = 1. 

( s \ 

The stationary volume V (X) is an increasing function of the compactivity X, since the 
powder "compressibility" 



dV is) (X) 1 

K ( x ) = — j^— = ^L 



dX 



Vi-v {s \x) 



(13) 



is proportional to the volume fluctuations over the ensemble of granular systems considered. 
Then, the compactivity X should be an increasing function of the vibration intensity T, 
because it is reasonable to expect that the powder became fluffier with higher vibration 
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intensities. This is, for instance, the behavior found in the simple two-volume model con- 
sidered in Ref . || . In general, this property must be checked once the relationship between 
T and X, X = f(T), has been derived for each specific model. 

(s) 

Of course, p\ is not a solution of the master equation when the intensity Y is time- 
dependent, and in general the system does not monotonically approach the steady distribu- 
tion. Nevertheless, define [compare with Eq. (0)] 



tf (s) (t)=£ft(*)m 



Pi(t) 



(14) 



where Pi(t) is again one solution of Eq. ([]]), and p\ s \x) depend on time through the com- 
pactivity X — X(t). If we define a statistical entropy as |T3| 



S*{t) = -X'£ Pi (t) ln Pi (t) 



(15) 



and use the notation 



V(t)=Y,ViPi(t) 



(16) 



for the actual average volume at time t, after a very simple algebra it is found that 



H (s \t) 



X 



V(t) - xs*(t)) - (v (s \x) - XS(X) 



(17) 



Thus H^ s \t) is proportional to the deviation of the actual "effective" volume ||,|5| at time t, 

Y(t) = V(t)-XS*(t), (18) 

from its stationary value. 

The time variation of H^ s ' is easily obtained as [12 



dH^ _ w _ v - Pi (t) dpf\x) dX 

{ ) r P P(x) dx dt> 



dt 



(19) 



where A^(t) is given by Eq. @, but replacing p'i(t) by p\ s \x). Taking into account Eqs. 
(10), it is found 



dt 



1 dX 

x 2 ~~dT 



V(t) - v {s \x) 



(20) 



Equation (EI]) does not have a well-defined sign and, therefore, the stationary distribution is 
not monotonically approached, in general, when the vibration intensity is time-dependent. 
However, in those processes such that the compactivity (or, equivalently, the vibration in- 
tensity r) increases monotonically in time, only the term 



B(t) 



1 dx v {s \x) 



X 2 dt 



(21) 
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is positive. Using the analogy between temperature and compactivity, we will refer to those 
processes as "heating" processes. If in a given "heating" process it is verified that 



lim£(t) = 0, (22) 

t — >oo 

we can conclude 

lim < . 23 

Since H^ s \t) is bounded below, the only possibility is in fact the equality, which implies 

lim A (s) (t) = (24) 

t — >oo 

and 

1 v 

$Zjrir v ® = - (25) 

From Eq. (]24|) and the irreducibility of the stochastic process, it follows that in the long 
time limit, 

lim # (t) = lim^ s) (X), (26) 

t — >oo t — >oo 

i. e., the system goes to the steady curve for long enough times, if the "heating" program 
verifies the condition given by Eq. (|22]) . 

The following picture emerges for the evolution of a granular system submitted to vibra- 
tions with increasing intensity: Starting from an arbitrary initial condition, in a first step 
the system tends to a behavior which is determined by the law of variation of the vibration 
intensity T, and the initial condition has been forgotten. This implies the existence of a 
special solution of the master equation, called the "normal" solution, that is approached by 
all the other solution after an initial transient period. Besides, if the system is "heated" 
slowly, in the sense that Eq. ( p2"D is verified, the normal solution tends afterwards to the 
stationary curve. This picture is similar to the one found in some models of structural 
glasses |mp|,P|. 



III. HILBERT'S METHOD AROUND THE STEADY CURVE 

In this Section we will use Hilbert's method to derive a quite general form of the normal 
solution of Eq. (]!]) near the steady state curve. We will focus on the class of models for 
granular media considered in the previous section, but the results can be directly extended 
to any system with a "canonical" distribution describing the stationary state. 

Let us consider Eq. (JTj) , rewritten in the form |14j 



^ = W(t) P (t) , (27) 

where p is a vector (column matrix) whose elements are the probabilities Pi(t) of the i-th 
state of the system at time t, and W is a square matrix with elements Wij given by 
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Wii(t) = W l3 (t) - Sa £ W kj {t) . (28) 

k 

If the transition rates are time independent and the detailed balance condition is verified, 
solving Eq. fl27|) is equivalent to obtain the solution of the eigenvalue problem 

W<p(q) = -\(q)<p(q) , (29) 

with X(q) > for all q. The eigenvectors tp(q) are completed with the stationary distribution 
p( s \ which is an eigenvector of W corresponding to the null eigenvalue, i. e., 

Wp {s) = . (30) 

Besides, if the Markov process is irreducible there is only one stationary distribution with 



all its components positive |f[4 |. The matrix W is hermitian with the following definition 



for the scalar product of any two vectors a and b, 

(M) = £^- (31) 

i Pi 

If the Markov process is irreducible and the detailed balance condition is fulfilled for 
all times, Eqs. (|29|-|3T|) remain valid for time dependent transition rates. Of course, the 
eigenvalues and eigenvectors will depend on time in general. The usual situation is that 
the transition rates depend on time through an externally controlled parameter like the 
temperature in a thermal system or the compactivity X in a granular medium. Thus, we 
will write sometimes in the following W(X), X(q,X), (p(q,X) and p( s \X). 

Hilbert's method consists in solving the master equation by means of the iterative process 



W(t)p(°\t) = 0, (32a) 

W(t)p (n) (t) = ^"fo® , n>l. (32b) 
In this way we obtain a probability distribution 

oo 

PH(t) = J2P (n) (t)i (33) 

n=0 

which is a solution of Eq. (|27|). The first term p^°'(t) gives the "stationary" distribution, 

P < i ) \t)= P ( i\x), (34) 

where X stands for the value of the compactivity at time t, X — X(t). The Hilbert solution 
Ph, constructed following the above rules, is a "normal" solution, because it only depends 
on the external law of variation of the transition rates, and it does not refer to any specific 
initial conditions. Nevertheless, the range of validity of Paif) is limited in general because of 



the divergence of Hilbert's expansion [17]. From a physical point of view, this divergence is 
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connected to the fact that we are expanding p(t), which may describe a situation arbitrarily 
far from the steady state, around the stationary solution p( s '(X). 

In "heating" processes we have shown in Sec. [IT] that a normal solution of the master 
equation exists and, under very general conditions, it tends to the stationary solution for very 
high compactivity. Thus in the high compactivity limit Hilbert's method can be useful, since 
it is an expansion around the steady state. We will restrict ourselves to the first correction 
pD(i), because the difference between the probability distributions Pu{t) and p^(X) is 
expected to be small in the high compactivity limit. The normalization of p^ s \X) implies 
that 



P (S) (X) 



dX 



dX 



(35) 



and, therefore, 



pV(t) = t(X) 



dpW(X)dX 
dX dt 



(36) 



where T{X) is the inverse operator of W(X) in the space orthogonal to the equilibrium 
distribution p( s \X). Using Eq. (§), we obtain 



d£\x) 



dX 



X 2 



v,-v (s \x) 



p 



(«) 



By introducing a function 



where cp(q,X) is the eigenvector defined in Eq. ( p9j) , we can rewrite Eq. 
form as 



dp [s \X) 
dX 



X 2 



(37) 
(38) 

in a vectorial 
(39) 



The above expression follows from the completeness of the eigenvectors and the property 



<p{q,x) 



dpWjxy 

dX 



X 2 



^9i{q,x) 



v t -v (s \x) 



1 

X 



(40) 



The term proportional to V (X) vanishes as a consequence of the orthogonality of the 
eigenvectors p^ s \X) and cp(q,X) for all q. Then, to first order in deviations from the 
equilibrium curve we have 



p H (t)=p^(X) 



1 dX 

x^m 



^X~ l (q : X)aq:XMq : X) 



(41) 



The notation used reflects the fact that the right hand side depends on time only through 
the compactivity X of the system. 
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From Eq. (filD it is possible to evaluate the average values of the physical quantities of 
the system over Hilbert's distribution in the first order approximation. For instance, the 
mean value of the volume is 



V H (t) = V (S \X) - ^^EA-^X^X) , (42) 



1 dX 



X 2 dt 



This expression can be written in a more transparent way by using the following identity 
for the powder "compressibility" defined in Eq. (pf), 

*P0 = 4E£ 2 (^), (43) 
q 

and the expression for the average relaxation time of the volume found in linear relaxation 
theory (see the appendix) 

T(x) = E g x-\ q ,x)e( q ,x) 



In this way, it is found 



V H (t) = V (s \x)-^(X)r(X). (45) 



It is important to note that the detailed balance condition is a key point in our derivation 
of Eq. ([45"D , while the existence of the normal solution does not need detailed balance to be 
satisfied, but it follows if the Markov process is irreducible in the long time limit. On the 
other hand, Eq. ([45]) can also be applied to estimate the departure from the stationary curve 
in "cooling" processes, if the system is initially prepared in the stationary state. In other 
words, Eq. (^) should be applicable to any situation near the stationary curve. Therefore, 
it allows to calculate the first order in the deviation of the volume from its steady value over 
the normal curve when the system asymptotically tends to the stationary curve, as it is the 
case of "heating" processes. The same remark can be made about Eq. ( (4~i~l) for Hilbert's 
distribution pn{t)- It contains all the information needed to evaluate one-time properties 
over the normal curve in situations close to the steady state. 



IV. THE MODEL 



In this Section we will briefly review a simple model for the vibrocompaction of a dry 
granular system which has been recently introduced [13], 18]- We consider a one-dimensional 



lattice with N sites. Each site can be either occupied by a particle or empty, i. e., occupied 
by a hole, with the restriction that there cannot be two nearest neighbor holes (such a 
configuration would be unstable). A variable m 8 is assigned to each site i, taking the value 
rrii — 1 if the site is empty, and mj = if there is a particle on it. Then, a configuration of the 
system is fully specified by giving the values of the set of variables m = {mi, m 2 , • • • , tun}. 

The dynamics of the system is defined as a Markov process, and formulated by means 
of a master equation for the probability p(m, t) of finding the system in the configuration 



m at time t 0] 



10 



^-p(m, t) = J2 [W(m\m')p(m, t) - W(m'\m)p(m, t)] , (46) 

m 

where W(m\m') is the transition rate from state m' to state m. The possible transitions 
in the system can be classified in three groups. Firstly, there are transitions conserving the 
number of particles, corresponding to purely diffusive processes. Their rates are given by 

W(010|100) = W(010|001) = -a , (47) 

with a being a constant. These transition rates must be understood as it is usually done, as 
the transition rates between states connected by the given rearrangement. Only the variables 
of the set of sites involved in the transition are indicated in the notation. Secondly, there 
are also transitions increasing the number of particles, with rates 

W(010|101) = ^a, (48a) 

W(001|101) = W(100|101) = . (48b) 
Finally, the transition rates for those processes decreasing the number of particles are 

W(01010|00100) = ^a 2 , (49a) 



W(01010|01000) = W(01010|00010) = ^a 2 . (49b) 
The transition rates in Eqs. (|47|)-(f49"|) define an effective dynamics for tapping processes 



taking place in a previous more general model presented in Ref. |fLq| . The parameter a 
characterizes the tapping process completely. Note that we have rescaled all the expressions 
of the transition rates in Ref. Jl3| in such a way that the time scale t is a measure of 
the number of taps. For a = no transition is possible in the system and, therefore, the 
parameter a measures the intensity of the vibration. The configuration with the highest 
density, i. e., no holes present, is completely isolated from the rest of the states, in the sense 
that no transition is possible from or towards it for any value of a. Aside from this particular 
state, all the other possible configurations are connected through a chain of transitions with 
non-zero probability for a ^ 0, and the Markov process defining the dynamics is then 
irreducible. 

When a does not depend on time, i. e., the intensity of vibration is constant, there is 
a unique stationary solution p^ s \m) of the master equation fl46|) . We will use the notation 
m O) f or a configuration of the system with k holes. Clearly, it is verified that 1 < k < 



(N + l)/2. It is easily found that (13 



-k/X 

p W( TO W) = « , (50) 
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where Z is the normalization constant, and X is related to the vibration intensity a through 
the relation 



a = e~ 1/x . (51) 

Comparison with Eq. (§) shows that in our model the number of holes k plays the role of the 
volume of the system and X is the compactivity of Edward's theory of powders . More 
precisely, k would be proportional to the excess volume from the densest state. Interestingly, 
the compactivity X is an increasing function of a, as it has been argued on general basis 
in Section |I|. Besides, the compactivity vanishes for the no vibrated case. The partition 
function can be analytically derived in the infinite system limit, with the result [|13[ 



InC = -^ lnZ = -hi2 + ln [l + (1 + 4ct) 1/2 ] . (52) 
From here, the stationary value of the density of holes is obtained in the standard way, 



x 



k" d\n( 1 

1 ~~N~ ~d(l/X) ~ 2 



l-(l + 4ay 1/2 . (53) 



(s) 

The stationary density of holes increases monotonically from the densest state x\ = to 
the fluffiest state x± = 1/2 as the vibration intensity increases monotonically from a = 
to a = oo, since the analogous to the "compressibility" 

dx[ s) e~V* 

K = —^r = ( 54 ) 

dX X 2 (l + 4e- 1 A) 3 / 2 

is positive definite for any value of X. Nevertheless, the densest configuration cannot be 
actually reached, since for a = no transition is possible and the system gets trapped in 
its initial state. As the number of holes k is an upper bounded variable, the compactivity 
X can take negative values, corresponding to configurations fluffier than those of positive 
compactivities. For instance, X — > 0~ gives a — > oo and Xi = 1/2, the least dense state. 

Finally, let us qualitatively study the dynamics of the system for time-independent tran- 
sition rates in the low vibration limit, a 1 or X — > + . In that limit, the stationary value 
of the density of holes x[ s ^ is very small, and a typical configuration of the system consist 
of a few holes, separated by long arrays of particles. Since x± ~ a, the mean distance 
between holes is of order a -1 . Moreover, at least in linear relaxation theory and at low 
compactivities, the evolution of the system will be mainly associated to diffusive processes. 
The characteristic time f of the relaxation process would be the square of the mean distance 
between holes divided by an effective diffusion coefficient, which can be estimated from Eq. 
([47|) as of the order a. Therefore, 

focoT 3 . (55) 

This result will be important in the following, since an estimate of the relaxation time in 
the linear relaxation regime is necessary when studying the time-dependent rates case, as it 
follows from Eq. 
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V. PROCESSES WITH TIME-DEPENDENT VIBRATION INTENSITY 



Next we will study the dynamical behavior of the model described in the previous section, 
when submitted to processes in which the vibration intensity a changes in time in a given 
way. Due to the relationship between a and the compactivity X, Eq. (|5TD, this is equivalent 
to consider that the compactivity X varies in time following a given law. As already men- 
tioned, we will refer to a process as a "heating" ( "cooling" ) one when the vibration intensity 
is monotonically increased (decreased). Such kind of processes have been already considered 
in the literature, both in real granular systems [|]|| and in simple models ||19|| . 

The general results of Sec. [Ij are applicable to this particular model. The only limitation 
is due to the loss of irreducibility of the dynamics for a = 0. Thus the existence of a 
special solution of the master equation, such that all the others approach it in the long 
time limit, applies to any program of variation of a except for "cooling" processes up to 
X = 0. As a consequence, for "heating" processes there will be a special "normal" curve, to 
which all the other solutions tend at a first stage. Later on, the system will approach the 
"stationary" distribution p^[a(t)] in the long time limit, provided that the condition in Eq. 
f(22T) is verified, i. e., 



Taking into account that is upper bounded by 1/2, and the relationship between a and 
X, Eq. (|5l]), the above condition can also be written as 

din aft) 

lim = . 57 

t^oo dt K ' 

Equation fJ57|) expresses a restriction for the "heating" programs driving the system to the 
stationary curve in the long time limit, but it does not affect at all the existence of the 
normal solution, which only depends on the ergodicity of the process. 
Application of Eq. (f45|) to the present model yields 



Xl (t) = x[ s) [X(t)} - d ^K[X{t)\r[X{t)\ + • • • , (58) 

where t(X) is the mean relaxation time of the density in the linear relaxation approximation. 
Upon writing the above expression we have taken into account that in our model the number 
of holes plays the role of the volume. For high compactivities the second term on the right 
hand side of Eq. (^) is negligible against the first one, since the compressibility k — > 
when X — > oo. Then, for high vibration intensities the system remains over the stationary 
curve. As the compactivity decreases, the system departs from that curve, as a consequence 
of the increase of the mean relaxation time r, which is expected to be proportional to the 
characteristic time f estimated in Eq. ([55]). 

There is some freedom when choosing the law of variation of the vibration intensity a. 
Our choice will be motivated by simplicity, but also by the analogies with a glass-like behavior 
previously found in granular systems PJTPfl. In supercooled liquids the temperature is usually 
varied at a constant rate |7]], so we have considered processes in which the compactivity 
changes linearly in time, 
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(59) 



with r > 0, which is equivalent to 

^ = ±ra(lna) 2 , (60) 

the plus sign corresponding to "heating" processes and the minus sign to "cooling" programs. 

The rest of this section is organized as follows. First, we study "cooling" processes. 
The system is initially put in the stationary state corresponding to a given value «o of the 
vibration intensity. Then, the compactivity X = — 1/ In a is decreased following Eq. fl59l). 
The existence of a phenomenon analogous to the laboratory glass transition of supercooled 
liquids arises. Afterwards, "heating" processes are considered, paying special attention to 
the appearance of hysteresis effects, and relating them to the trend of the system to approach 
the normal curve. 



A. "Cooling" processes 

We consider the continuous decreasing of the compactivity of the system from a given 
initial value X Q down to X = 0. The latter corresponds to a = 0, i. e., no vibration. The 
system is initially placed in the stationary state corresponding to the value a = exp(— 1/X ) 
of the vibration intensity. Then, the compactivity is decreased following the law 

dX 

-dT = ~ rc - (61) 



Our starting point will be Eq. Q58D, particularized for the "cooling" process we are 
considering, i. e., 

Xl (t) = x?[X(t)} + r c K [X(t)] r[X(t)} + 0(r c 2 ) . (62) 

As we have already discussed, from the above equation follows that the system remains in 
equilibrium at "high" compactivities, since the second term in the right hand side vanishes 
for X — > oo. Nevertheless, as the compactivity becomes smaller this term grows, due to the 
increase of the relaxation time r. This means that there will exist a range of values of the 
compactivity in which the second term is comparable to the first one. A rough estimate 
of the value of the compactivity X g at which the system would depart from the stationary 
curve can be obtained by equalling both terms. If we consider that the system is slowly 
"cooled", r c <C 1, it is also a g = e~ l l Xg <C 1, and we can approximate both terms for their 
leading behaviors in that limit, namely 

x[ s) {X)~a } (63a) 



K (X) = ~ a (In of (63b) 
dX 



dx[ s ) 

1x 

t(X) ~ r a- :i (63c) 
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where tq is a constant of the order of unity. Therefore, we get 



r c r (In a g ) = a g . (64) 

For r c <C 1 it is 

f 3 |lnr c | 2/3 . (65) 



OCg 



We have omitted the factors containing r , as well as any other factor of the order of unity. 

For a < a g the system is effectively "frozen" , due to the divergent tendency of the relax- 
ation time no more transitions are possible, and the density of holes will be approximately 
constant in this region. A measure of the effective number of transitions left to the system 
before reaching a = from a given time t is given by the scale |T0|JT6[| 



s(t) = C dt 'WW)y (66) 



where t is the time instant for which the "cooling" program finishes, i.e. a(t ) = 0. In 
that way, the system would get frozen for t > tf such that s(tf) = 1, being easily obtained 
that a(tj) ~ a g . Using Eq. (|65"D it is possible to estimate the leading order value of the 
compactivity at which the system gets frozen, 

^ 9 \na g | In r c | ' ^ 
This kind of behavior has been also found numerically in a granular system model |19 



The inverse logarithmic dependence of the cooling rate is typical for the laboratory glass 
transition temperature of supercooled liquids f7j, and it has also been analytically derived 
in some simple models of structural glasses pT0|JT6[ . 



Since the density of holes remains nearly constant for a < a g , and the two first terms of 
Hilbert's expansion ( |5"2"D are of the same order for a ~ a g , it is reasonable to expect that 

Xijes = Km Xi(t) oc x < f\a g ) = a g ~ r c 1/3 |lnr c | 2/3 , (68) 

t— >to 

where x^res is the residual value of the density of holes, extending again to this model of 
granular system the terminology of structural glasses. 

In order to check the above results, Fig. |l] shows the residual value of the density of holes 
as a function of the "cooling" rate, measured by the parameter 

5 = r c (\nr c ) 2 , (69) 



in a log-log scale. The numerical result agrees with the theoretical prediction, Eq. 
since the curve is well fitted by a straight line with a slope approximately equal to 1/3. 
In Fig. Q the evolution of the density of particles, p — 1 — xi, in a "cooling" process with 
rate r c = 10~ 5 is plotted. For comparison, the equilibrium curve, given by Eq. (|53|) , is also 
shown. The estimate of the freezing compactivity from Eq. ( |6TD is X g ~ 0.26, which is seen 
to be in good agreement with the region in which the Monte Carlo density is approximately 
constant. Similar behaviors have been observed for other small values of the cooling rate. 
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B. "Heating" processes and hysteresis effects 

Let us analyze "heating" processes, i. e., processes in which the vibration intensity is 
monotonically increased. In this kind of processes, the dynamics of the system is irreducible. 
Thus there is a "normal" solution of the master equation, such that any other solution tends 
to it. Besides, if the heating program verifies Eq. (|57|), the normal solution approaches the 
stationary curve for large enough times. We are going to discuss how these results can be 
applied to our model, in order to understand its behavior when "heated" from a = 0. 

The compactivity will be increased according to the law 

dX . 

where is the rate for this process. According to Eq. (plf), the system will approach the 
stationary curve in the high compactivity (long time) limit. In the vicinity of the stationary 
curve, the evolution of the system is given by 

Xl (t) = x[ s) [X(t)} - r h K[X(t)}r[X(t)} + Otf) . (71) 

This equation explains why, in "heating processes" , the system tends to the stationary curve 
following a curve different from the "cooling" one, even when r c = r^. It follows directly 
from the comparison of Eqs. ( |7T|) and Eq. (0), by noting that the deviation from the 
stationary behavior is of opposite signs in "heating" and "cooling" processes. Therefore, 
hysteresis effects show up. 

However, perhaps the main result for "heating" processes is the existence of the normal 
solution. Fig. |] shows the evolution of the density of particles in a heating process with 
rh = 10~ 5 . The shape of the curve depends on the initial condition for X = 0. Two different 
initial preparations of the system have been considered. The system was previously cooled 
down to X = 0, following two different linear programs with r c = 10~ 5 and r c = 10~ 3 , 
respectively. For the sake of clarity, these "cooling" curves are not shown. From Fig. |3] it is 
seen that both "heating" curves tend to a common behavior and, afterwards, they approach 
the stationary curve for high compactivities. Also plotted is the normal curve, which was 



obtained by starting from the loosest packing state, x\ = 0.5 |20 . 

Figure [| depicts a particular cycle of "cooling" and "heating" with the same rate, namely 
r c = rh = 10 -5 , as well as the normal curve of the "heating" process. A behavior similar 



to the one found in real granular systems and also in the "Tetris" model is 

observed. When starting the heating process from the loosest packing state the normal 
curve is obtained, which tends to the stationary behavior in the high vibration intensity 
limit. Afterwards, cooling and reheating with the same rate leads to the other two curves 
of the figure. These are approximately "reversible" for very small rates, since the deviation 
from the stationary curve is smaller the smaller the rate. Nevertheless, they cannot be used 
to obtain the stationary values of the density for low compactivities, due to the glass-like 
kinetic transition. On the other hand, at high compactivities the deviations of the system 
from the stationary curve for the "cooling" and the "heating" processes are symmetric, as 
predicted by Eqs. (|62|) and (fH|). 
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VI. FINAL REMARKS 



In the first part of the paper we have considered a wide class of models of granular 
systems, namely those models whose dynamics under tapping is described by means of a 
master equation. It has been assumed that the tapping process is such that the system is able 
to explore the whole space of metastable states, when it is tapped with any vibration intensity 
r ^ 0. Then, if T changes in time, and the Markov process remains irreducible in the long 
time limit, all the solutions of the master equation tend to approach a special solution, called 
the "normal" solution for the given vibration program. Besides, if the stationary distribution 
for constant T is consistent with Edward's statistical mechanics theory of powders, the 
normal solution tends to the stationary curve in the high vibration intensity regime, provided 
that the "heating" program is not too fast, in the sense given by Eq. (T2). It has also been 
argued that the compactivity X of Edward's theory should be an increasing function of the 
vibration intensity V . A quite general prediction for the behavior of the normal solution 
in situations near the steady state has been found, by using Hilbert's method to solve the 
master equation. This result can be applied to any system having a "canonical" steady state 
distribution. 

Very recently, a simple model for a granular system submitted to vertical vibration with 
an stationary state described by Edward's theory of powders has been introduced [|i"3| , |TB|| . 
The model belongs to the class of systems described in the previous paragraph, and here 
we have investigated its behavior in processes in which the vibration intensity depends on 
time. Such kind of processes have been carried out in real granular systems plfl, and also 



investigated in some models |19| . For the sake of concreteness we have taken the compactivity 
as a linear function of time. 

The behavior of the model under "cooling" processes up to zero vibration intensity 
exhibits a phenomenon similar to the laboratory glass transition. The "cooling" evolution 
departs from the stationary curve, and freezes, in a narrow region around the value of 
the compactivity X g such that the effective number of transitions left to the system until 
reaching X = becomes smaller than unity. This allows us to estimate the dependence of 
X„ and the residual values of the density upon the "cooling" rate. The results are analogous 



to those obtained previously in simple models of structural glasses [10,16 



In "heating" processes, a crucial role is played by the "normal" solution of the master 
equation, which is completely determined by the program of variation of the compactivity, 
and attracts any other solution, independently of the initial condition. The hysteresis effects 
found when the system is "cooled" and "reheated" are related to the trend of the system 
to approach the normal curve. The normal curve corresponds to "heating" the system from 
the loosest state. In a first stage, the system tends to approach monotonically the "normal" 
curve and, for longer times, corresponding to high enough compactivities, the stationary 
curve is reached if the "heating" program verifies Eq. (p>7j). 

The work reported here suggests a very close relationship between structural glasses and 
the behavior of a granular system under vibrations, supporting previous results in the same 
direction [|1^,|1^,^] . It also raises some questions, both from a theoretical and experimental 
point of view. We think that it is worth checking the existence of the normal solution in real 
granular systems, and whether it can be constructed starting from the loosest packing state. 
Also, it seems interesting to know if the glass analogy also extends to the linear relaxation 
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regime, giving a stretched exponential behavior of the response functions, or if the inverse 
logarithmic relaxation law @ remains valid in such a region. The latter possibility would 
lead to the need of investigating the reason for such a behavior. 
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APPENDIX A: AVERAGE RELAXATION TIME OF THE VOLUME 

Linear relaxation describes the evolution of the system at constant value of the com- 
pactivity X, when starting from the steady state corresponding to a compactivity X + AX, 
with AX very small, in the sense that we can approximate 

p {s \X + AX) = pW(X) + ^^ AX . (A1 ) 

aX 

Taking into account Eq. (0), the evolution of the probability distribution from this initial 
condition, p(t = 0) = p^(X + AX), is 



A X 

pit) = p^(X) + —Y,^X)^X)e- tx ^) , (A2) 



q 



because ip(q, X) is the eigenvector of the evolution operator corresponding to the eigenvalue 
X(q,X). From the above expression it is straightforward to calculate the time evolution of 
the average volume, 

V(t) -F (S) (X) = §Ee 2 (a)e- (AM , (A3) 

i 

where we have made use of the definition of the function £(q,X), Eq. (j3§|). 

The volume linear relaxation function corresponding to a value of the compactivity X is 
defined in the usual way, 

(A4, 

V(t = 0)-V [ >(X) 

yielding 

Mt ' X)= £,?(<,, x) ■ (A5) 

On the other hand, the linear relaxation time for the volume r(X) is given by definition by 
the area under the curve <pv(t] X) as a function of t. This leads to Eq. (H). 



18 



REFERENCES 



[1] J. B. Knight, C. G. Frandich, C. N. Lau, H. M. Jaeger, and S. R. Nagel, Phys. Rev. E, 
51, 3957 (1995). 

[2] E. R. Nowak, J. B. Knight, E. Ben-Naim. H. M. Jaeger, and S. R. Nagel, Phys. Rev. E 
57, 1971 (1998). 

[3] H. M. Jaeger, Proceedings of the NATO/ASI Workshop on Dry Granular Materials, 

Cargese 1997 (unpublished). 
[4] S. F. Edwards and R. B. S. Oakeshott, Physica A 157, 1080 (1989). 
[5] A. Mehta and S. F. Edwards, Physica A 157, 1091 (1989). 
[6] S. F. Edwards and D. V. Grinev, Phys. Rev. E 58, 4758 (1998). 
[7] G. W. Scherer, Relaxation in Glass and Composites (Wiley, New York, 1986). 
[8] M. J. Ruiz-Montero, Modelos sencillos de relajacion estructural en vidrios, PhD Thesis, 

Universidad de Sevilla, 1992. 
[9] S. Franz and F. Ritort, Europhys. Lett. 31, 507 (1995). 
[10] J. J. Brey and A. Prados, Phys. Rev. E 49, 984 (1994). 

[11] A. Prados, Fenomenos de relajacion en modelos sencillos que muestran comportamiento 
vitreo, PhD Thesis, Universidad de Sevilla, 1993. 

[12] J. J. Brey and A. Prados, Phys. Rev. E 47, 1541 (1993). 

[13] J. J. Brey, A. Prados and B. Sanchez- Rey, Physica A 275, 310 (2000). 

[14] N.G. van Kampen, Stochastic Processes in Physics and Chemistry (North-Holland, Am- 
sterdam, 1992). 

[15] J. J. Brey and A. Prados, Phys. Rev. A 42, 765 (1990). 

[16] A. Prados, J. J. Brey and B. Sanchez-Rey, Phys. Rev. B 55, 6343 (1997). 

[17] P. Resibois and M. de Leener, Classical Kinetic Theory of Fluids (Wiley, New York, 
1977). 

[18] J. J. Brey, A. Prados and B. Sanchez-Rey, Phys. Rev. E 60, 5685 (1999). 
[19] M. Nicodemi and A. Coniglio, Phys. Rev. Lett. 82, 916 (1999). 

[20] It can be easily understood why the usual choice for the initial condition made in 
models of structural glasses |T^]T2]]T6[[ does not work in the present context. In those 
models the initial condition was taken to be the limit of the stationary curve for zero 
temperature, which corresponds here to the densest state, x± = 0. But the densest 
configuration does not evolve for any value of the compactivity, as it should also happen 
in a perfectly ordered real granular system. Then, in our model the normal solution 
corresponds to putting the system initially in the other extreme limit of the density, the 
loosest configuration. 

[21] L. P. Kadanoff, Rev. Mod. Phys. 71, 435 (1999). 



19 



FIGURES 



FIG. 1. Residual value £i,res of the hole density as a function of the "cooling" parameter S, 
obtained from Monte Carlo simulation of the system. The continuous line is the best linear fit, 
having a slope equal to 0.32, very close to the theoretical value. 

FIG. 2. Evolution of the particle density p with the compactivity X for a cooling process 
with rate r c = 10 (diamonds). The dotted line is the equilibrium function. The estimate for the 
freezing compactivity X g is also shown. 

FIG. 3. Evolution of the density of particles p with the compactivity X for a heating process 
with rate = 10~ 5 . Starting from the top, the symbols correspond to heating processes after 
coolings with rates r c = 10 -3 (crosses) and r c = 10~ 5 (diamonds). Both curves tend to approach 
the normal curve (continuous line), whose initial condition is the loosest configuration, p = 0.5. 

FIG. 4. Reversible-irreversible cycle, corresponding to a rate r c = = 10 -5 . Starting from 
the loosest configuration, the system is heated (squares), tending to the equilibrium curve (dotted 
line) for high enough compactivities. Afterwards the system is cooled (diamonds) and reheated 
(crosses). The last two processes are approximately reversible. 
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